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Abstract 

We develop R and C code for Individual Difference Multidimensional Scaling, both 
for the INDSCAL and the IDIOSCAL case. In addition to the new SMACOF algorithms 
to minimize the stress loss function we use expressions for the second derivatives with 
respect to the group configuration and the individual weights to compute perturbation 
ellipsoids. 
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version, the bib files, the complete Rmd hie with the code chunks, and the R and C source 
code. 


1 Introduction 

In Individual Differences Multidimensional Scaling (IDMDS) we minimize the stress loss 
function, first proposed by Kruskal (1964a), Kruskal (1964b), and defined as 

-1 m 

<r(x u ■ • •, x m ) = - E EE w ijk (s ijk - d l3 {x k )f (i) 

2 k =11 <i<j<n 

over the n x p configurations X\, ■ ■ ■ , X rn . In (1) the Wij k are hxed and known weights, the 
Sijk are hxed and known dissimilarities, and d t3 (X k ) is the Euclidean distance between rows 

1 and j of X k . 

The conhgurations are constrained to be of the forrriX fc = XT k , where T k is a square matrix 
of order p. In IDMDS we typically choose T k to be either a general matrix (IDIOSCAL) or a 
diagonal matrix (INDSCAL). 

2 Weighted Euclidean Distances 

The Euclidean distance d t] (X k ) can be expressed in various ways. Some interesting ones, for 
our purposes, are 


dij(X k ) y(x* rj ) C k (xi Xj ) 

= x ' k A ij X k 
= \Jx'(C k ® Aij)x, 

where 

C k = T k T' k , 

Aij = (Cj ~ e j)( e i — e j) > 

Xi = X'ei, 
x = vec(X). 

In these expressions the e t and e 3 are unit n-element vectors (columns of the identity matrix), 
with a single element equal to one and all other elements equal to zero. Thus Ai 3 is an n x n 
matrix, with elements (i,i) and (j, j) equal to +1, elements ( i,j ) and j,i ) equal to —1, and 


( 2 ) 

( 3 ) 

( 4 ) 
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all other elements equal to zero. We use (g) for the Kronecker product. The matrix C k (g) 
is a symmetric rip x rip matrix, built with the p x p blocks CkstAij of size n x n. 

Instead of expressing the weighted distance as the square root of a quadratic form in x, we 
can also find a similar representation as the square root of a quadratic form in the T k . 

dpx t ) = ytr TppA~XT k = v'iyt® XV)‘* (5) 

Here t k = vec(T*.). The matrix / (g) X'A VJ X is the direct sum of p copies of the p x p matrix 
X'AijX = (xi — Xj){xi — Xj)'. If the T k are diagonal, as in INDSCAL, then we have a more 
efficient representation. 


dij(X t ) = ytr nX'A tJ XT t = ^X'^Xt,, (6) 

where t k is now the diagonal of T k . 

3 SMACOF 

The basic theory for SMACOF algorithms is in De Leeuw (1977). Applications to IDMDS 
are in De Leeuw and Heiser (1977) and De Leeuw and Heiser (1980). The IDMDS algorithm 
suggested in these papers, which is implemented in De Leeuw and Mair (2009), is different 
from the one we propose in this section. 

Let’s start by supposing, without loss of generality, that the dissimilarities are normalized, in 
the sense that 


-i m 

aEEE Wijktfjk = L ( 7 ) 

k= 1 l<i< 7 <n 

The IDMDS loss function becomes 

m -t m 

<r(x u • • •, x m ) = 1 -EEE w ljk s ijk d tJ (x k ) + -EEE Wi jk d%{x k ). (s) 

k =1 k =1 l<i<7<n 


Now 

(9) 
( 10 ) 

and thus 


ddij(X k ) 

dx 

ddUX k ) 


dx 


d’ij (A k 
= 2 (C k (g) Aij)x, 


(C k (g) A i: j)x, 


3 



^( 4 - ,x m ) 

dx 


V+x-B^x, 


where 


K =E( C k® V *)> 

k= 1 

m 

B*{x) = E( C k ® B k (x)), 

k= 1 

and 


(ii) 


( 12 ) 

(13) 


EE (14) 

Bk{x) =EE w vk ~7 7Y~ , Aij. (15) 

The SMACOF iterations to minimize (1) over x for given T k are 

= V+B^x^xW, (16) 

where v is the iteration counter and superscript + is the Moore-Penrose inverse. 

Now we do the same for the T k . First the 1D10SCAL case 


ddjji X k ) 
dt k 

dd%{X k ) 

dt k 


d ij (X t ) i ' I ® X ' AilX)tk ' 
2(1 ® X'Atj X)t k , 


( 17 ) 

(18) 


and thus 


° a(Xl \''' ,Xm) = (I ® X'V k X)t k - (./® X'B k {x)X)t k . 

dtk 

The SMACOF iteration step to update T k for fixed X is 

4" +1) = (I <D (X'V t X)+)(I ® X'B t (x)X)t£K 


(19) 


( 20 ) 
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Note that this means we can update each column of T k separately while, of course, we were 
already updating each T k separately. 

For INDSCAL we have an even simpler iteration on the diagonal of T k 

4" +1) = {X'V h X)+X'B k {x)Xt%\ (21) 

The overall SMACOF algorithm, implemented in the smacof IDI0SCALO and smacof INDSCAL () 
functions in the appendix, alternates iterations to improve X for fixed T k with iterations to 
improve the T k for fixed X. In our actual implementation we only use a single one of each of 
these two inner iterations to define a SMACOF step, and after each of these SMACOF steps 
we test for convergence. Of course the general convergence theory supports algorithms with 
more than one inner iteration step, but we have not experimented with varying that aspect 
of the algorithms. 


4 Example (Rectangles) 

Our example uses data from Borg and Leutner (1983). The data are included in the smacof 
package. In this example m — 2, p — 2, and n = 16. The help hie says 

42 subjects are assigned to two groups of 21 persons. 120 stimulus pairs of rect¬ 
angles are presented. For the first group (width-height; WH), the rectangles were 
constructed according to a design as given in rect constr. For the second group 
(size-shape; SS) the rectangles were constructed according to a grid design, which 
is orthogonal in the dimensional system reflecting area (size), and width/height 
(shape). All subjects had to judge the similarity of the rectangles on a scale from 
0 to 9. 

The IDIOSCAL and INDSCAL iterations are both started using the output of the idioscal () 
and indscalO functions in the smacof package as initial configurations. For IDIOSCAL 

h.idio <- smacofIDIOSCAL (w, delta, res.idio$cweights, res.idio$gspace) 


## 

itel 

1 

sold 

Inf 

sa 

197.115996 

sb 

0.070464 

## 

itel 

2 

sold 

0.070464 

sa 

0.061518 

sb 

0.060751 

## 

itel 

3 

sold 

0.060751 

sa 

0.057689 

sb 

0.057362 

## 

itel 

4 

sold 

0.057362 

sa 

0.056299 

sb 

0.056118 

## 

itel 

5 

sold 

0.056118 

sa 

0.055761 

sb 

0.055650 

## 

itel 

6 

sold 

0.055650 

sa 

0.055531 

sb 

0.055458 

## 

itel 

7 

sold 

0.055458 

sa 

0.055417 

sb 

0.055369 

## 

itel 

8 

sold 

0.055369 

sa 

0.055355 

sb 

0.055322 

## 

itel 

9 

sold 

0.055322 

sa 

0.055317 

sb 

0.055294 

## 

itel 

10 

sold 

0.055294 

sa 

0.055292 

sb 

0.055277 

## 

itel 

11 

sold 

0.055277 

sa 

0.055276 

sb 

0.055266 

## 

itel 

12 

sold 

0.055266 

sa 

0.055266 

sb 

0.055259 

## 

itel 

13 

sold 

0.055259 

sa 

0.055258 

sb 

0.055254 
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## itel 

14 sold 

0.055254 sa 

0.055254 sb 

0.055250 

## itel 

15 sold 

0.055250 sa 

0.055250 sb 

0.055248 

## itel 

16 sold 

0.055248 sa 

0.055248 sb 

0.055247 

## itel 

17 sold 

0.055247 sa 

0.055247 sb 

0.055246 

## itel 

18 sold 

0.055246 sa 

0.055246 sb 

0.055245 


For INDSCAL the iterations are 

h.diag <- smacofINDSCAL (w, delta, res.diag$cweights, res.diag$gspace) 


## 

itel 

1 

sold 

Inf 

sa 

197.116619 

sb 

0.070504 

## 

itel 

2 

sold 

0.070504 

sa 

0.062663 

sb 

0.061343 

## 

itel 

3 

sold 

0.061343 

sa 

0.058372 

sb 

0.057804 

## 

itel 

4 

sold 

0.057804 

sa 

0.056646 

sb 

0.056361 

## 

itel 

5 

sold 

0.056361 

sa 

0.055928 

sb 

0.055770 

## 

itel 

6 

sold 

0.055770 

sa 

0.055611 

sb 

0.055517 

## 

itel 

7 

sold 

0.055517 

sa 

0.055459 

sb 

0.055399 

## 

itel 

8 

sold 

0.055399 

sa 

0.055378 

sb 

0.055339 

## 

itel 

9 

sold 

0.055339 

sa 

0.055331 

sb 

0.055305 

## 

itel 

10 

sold 

0.055305 

sa 

0.055302 

sb 

0.055284 

## 

itel 

11 

sold 

0.055284 

sa 

0.055283 

sb 

0.055271 

## 

itel 

12 

sold 

0.055271 

sa 

0.055271 

sb 

0.055262 

## 

itel 

13 

sold 

0.055262 

sa 

0.055262 

sb 

0.055257 

## 

itel 

14 

sold 

0.055257 

sa 

0.055257 

sb 

0.055253 

## 

itel 

15 

sold 

0.055253 

sa 

0.055253 

sb 

0.055250 

## 

itel 

16 

sold 

0.055250 

sa 

0.055250 

sb 

0.055248 

## 

itel 

17 

sold 

0.055248 

sa 

0.055248 

sb 

0.055247 

## 

itel 

18 

sold 

0.055247 

sa 

0.055247 

sb 

0.055246 


As the minimum of the stress already suggests the two solutions for A" are basically indistin¬ 
guishable. We plot the one from INDSCAL. 
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INDSCAL Rectangles 



diml 

Since there are only two different distance matrices we do not plot the T} ; . For IDIOSCAL 
they are 

## [[ 1 ]] 

## D1 D2 

## D1 0.9969528775 -0.1401381035 
## D2 -0.1900447368 1.0747955535 

## 

## [[ 2 ]] 

## D1 D2 

## D1 1.0810614200 0.1374772239 
## D2 0.1829126197 0.8261765584 

and for INDSCAL we have 

## [[ 1 ]] 

## [, 1 ] [, 2 ] 

## [1,] 0.8760383307 0.000000000 
## [2,] 0.0000000000 1.312868518 
## 

## [[ 2 ]] 

## [, 1 ] [, 2 ] 

## [1,] 1.162505002 0.0000000000 
## [2,] 0.000000000 0.8206214109 

The Tfc from IDIOSCAL are hevily diagonally dominant, which explains why INDSCAL and 
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IDIOSCAL are very similar. 


5 Example (Colors) 

Our second, somewhat larger, example uses data from Helm (1959), also included in the 
smacof package. From the help file: 

A detailed description of the experiment can be found in Borg and Groenen 
(2005), p. 451 with the corresponding Table 21.1. containing distance estimates 
for color pairs. There were 14 subjects that rated the similarity of colors, 2 of 
whom replicated the experiment. 10 subjects have a normal color vision (labelled 
by N1 to N10 in our list object), 4 of them are red-green deficient in varying 
degrees. In this dataset we give the dissimilarity matrices for each of the subjects, 
including the replications. 

The IDIOSCAL iterations are: 

h.idio <- smacofIDIOSCAL (weights, delta, idi,helm$cweights, idi,helm$gspace) 


## 

itel 

1 

sold 

Inf 

sa 

632.304054 

sb 

0.063500 

## 

itel 

2 

sold 

0.063500 

sa 

0.043978 

sb 

0.043829 

## 

itel 

3 

sold 

0.043829 

sa 

0.043316 

sb 

0.043283 

## 

itel 

4 

sold 

0.043283 

sa 

0.043092 

sb 

0.043083 

## 

itel 

5 

sold 

0.043083 

sa 

0.043015 

sb 

0.043013 

## 

itel 

6 

sold 

0.043013 

sa 

0.042989 

sb 

0.042988 

## 

itel 

7 

sold 

0.042988 

sa 

0.042980 

sb 

0.042980 

## 

itel 

8 

sold 

0.042980 

sa 

0.042977 

sb 

0.042977 

## 

itel 

9 

sold 

0.042977 

sa 

0.042976 

sb 

0.042976 


The INDSCAL iterations are: 

h.diag <- smacofINDSCAL (weights, delta, ind.helm$cweights, ind.helm$gspace) 


## 

itel 

1 

sold 

Inf 

sa 

631.986429 

sb 

0.065305 

## 

itel 

2 

sold 

0.065305 

sa 

0.046937 

sb 

0.046808 

## 

itel 

3 

sold 

0.046808 

sa 

0.046355 

sb 

0.046325 

## 

itel 

4 

sold 

0.046325 

sa 

0.046131 

sb 

0.046123 

## 

itel 

5 

sold 

0.046123 

sa 

0.046037 

sb 

0.046035 

## 

itel 

6 

sold 

0.046035 

sa 

0.045995 

sb 

0.045993 

## 

itel 

7 

sold 

0.045993 

sa 

0.045974 

sb 

0.045973 

## 

itel 

8 

sold 

0.045973 

sa 

0.045963 

sb 

0.045962 

## 

itel 

9 

sold 

0.045962 

sa 

0.045956 

sb 

0.045956 

## 

itel 

10 

sold 

0.045956 

sa 

0.045953 

sb 

0.045952 

## 

itel 

11 

sold 

0.045952 

sa 

0.045951 

sb 

0.045950 

## 

itel 

12 

sold 

0.045950 

sa 

0.045949 

sb 

0.045949 

## 

itel 

13 

sold 

0.045949 

sa 

0.045948 

sb 

0.045948 


INDSCAL Helm Configuration 
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6 Stability 

In this paper we study, in addition, the stability of IDMDS solutions. We define stability 
in the usual way. In inferential statistics we look at the second derivatives of the likelihood 
function at the maximum likelihood estimate. In MDS we look at the second derivatives of 
stress at a local minimum. Our results extend the results in De Leenw (2017) from MDS to 
IDMDS. 

In general, if a function / : M n M is two times continuously differentiable then 

f(x) = f(y ) + {x - y)'Vf(y) + * (x - y)'V 2 f(y)(x - y) + o(||a? - y\\). 

If / has a local minimum in y then T>f(y) = 0 and V 2 f(y) > 0. Consider all x of the form 

x = Nj , (22) 

[V2 \ 

where y -2 has the last n — m elements of y. Then 

f(x) = f(y ) + *(zi - y 1 )’H(y){x l - y x ) + o(||xi - 2 /i||), (23) 

where H(y) is the leading m x m submatrix of the T> 2 f(y). Thus the set of all x of the form 
(22) with f(x) < (1 + e)f(y) can be approximated by choosing X\ in the ellipsoid 


Od - yi) , H(y)(x 1 - 7/i) < 2 ef(y) (24) 

We could choose e to be 0.1, for example, which means we look at a perturbation region 
where stress is at most 10% larger than the minimum we have found. 

To apply these general results to IDMDS we need the second derivatives of stress. Of course 
these second derivatives can also be used for other purposes, such as checking the necessary 
conditions for a local minimum, or for implementations of Newton’s method. Start with the 
derivatives with repect to X. We have 


d 2 d i:j (X k ) 1 j (C k ® Aij)xx'(Ck <S> Aij) 


and, of course, 


d 2 dj,(X k ) 

dxdx 


Ck & -^-ij • 


(25) 


(26) 


Thus 
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where 


dxdx 




( 27 ) 




= EEE “>y* 

k =1 


Oijk 

dij (X k 


(C k ® A 




v) 


Aij)xx'{C k 

w 


A, 




(28) 


We now look at the second derivatives with repect to the T fe . First the IDIOSCAL case, 
which has 


dt k dt 


I <g> X'AijX, 


and 


(29) 


9%(AV) 1 y , (I®X’A ll X)tMI®X’A,,X) 

am* <mvo V u ' 4(vo 

,Am) = (/ ® X'VtX) - G t (4), 

OtkOtk 

where 


(30) 

(31) 


G k (t k ) 


EE ^ijk 

1 


dijk 

dij (X k ) 


(/ ® XMy X) 


(/ ® X'AijXy^I ® X'AijX ) 
<%(**) 


(32) 


For the INDSCAL case we have the same formulas, but without the /(g) part. Thus 


where 


d 2 q(AV- - ,X m ) 
dt k dt k 


X'V k X-G k (t k ), 


(33) 


^(4) = E E w ^ k /(Yj f vM " v - • (34) 

We should perhaps mention that the perturbation results for IDIOSCAL are limited by the 
fact that the representation of X k as XT k is far from unique. Besides the obvious, and rather 
harmless, unidentifiability due to translation and expansion, we also have XT k = (A -Q)(Q~ 1 T k ) 
for any nonsingular Q. This suggests we should incorporate an identification condition such 
as T\ — I in our equations. We haven’t done this yet. Until we have chosen the identification 
condition the stability analysis for IDIOSCAL is not yet available. Note that the identification 
condition is not needed for the SMACOF algorithm, but it is needed for the stability analysis. 
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7 Example (Colours Stability) 


We now apply our second derivative information to the Helm INDSCAL analysis. First for 
the configuration X. 

par(pty="s") 

hh<-smacofDerivativesX (weights, delta, h.diag$b, h.diag$x) 
smacofEllipsesX (h.diag$x, hh$h, hh$s, .05, 1, 2) 


CM 

E 

T3 


CO 

O 

d 


o 

d 


o 

d 

i 


CO 

o 

d 

i 



-0.04 -0.02 0.00 0.02 0.04 


dim 1 

For the weights T k in INDSCAL we compute a separate 2x2 matrix of derivatives for each 
of the 16 replications. 

bbb <- t(sapply (h.diag$b, diag)) 
hh <- array (0, c(2, 2, 16)) 
for (k in 1:16) { 

hk<-smacofDerivativesTKDiag (weights[[k]], delta[[k]], h.diag$b[[k]], h.diag$x) 
hh[, , k] <- hk$h 

} 

smacofEllipsesTDiag (bbb, hh, h.diag$s, .002, 1, 2) 
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T3 


C\J 


E 



00 

o 



0.8 0.9 1.0 1.1 1.2 1.3 


dim 1 


It is remarkable that both the orientation and size of the ellipses are basically the same over 
the 16 replications. It remains to be seen if this is specific to this example or if it happens 
more generally. 

8 Projection 

In our stability analysis for X we change row x i} and we keep all other rows and all T\ fixed 
and constant. As an alternative we could look at 


o,{X) = min cr(Xi, ■ ■ • , X k ). 


(35) 


We project out the T *. and loss becomes a function of X only. Now 


d 2 a.(X) _ d 2 a( AV" ,X k ) 


dxdx dxdx 



We have not implemented this alternative defintion of stability yet, mainly because it requires 
in addition the mixed second partials with respect to x and the t k . Clearly at a minimum 


av.( x) < av(. n,--- ,x t ) 


( 37 ) 


dxdx 


dxdx 


which implies the perturbation ellipsoids will be larger. 
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9 Nonmetric MDS 


In our results so far we have treated the Sijk as fixed constants. In Nonmetric Multidimensional 
Scaling (NMDS), however, stress can be defined as 

^ m 

a(x h ---,x min EE u^ijk ( K ^ijk dtj (Ah) ) , (38) 

Z k=l 5keAk l<i<j<n 

which makes the dissimilarities a function of the Ah- This projection changes the perturbation 
results for the X in the same way as the result in the previous section did. But with NMDS 
there is an additional complication. Typically the dissimilarities are minimized over the 
polyhedral cones that define monotone regression. Although the partial minimum of stress 
over the cone is differentiable with respect to the X* ; , it is not twice differentiable, because 
the monotone regression transformation itself is not differentiable. Again, we have chosen 
not to implement the resulting complications (yet). 

The perturbation results we have can be applied in NMDS, but the perturbation regions 
must be interpreted in terms of moving the points in the configuration for a given set of 
dissimilarties, which can of course be the optimal dissimilarities computed by NMDS. 

10 Appendix: Code 
10.1 smacoflDIOSCAL.R 

# w is list of length m with n x n matrices 

# delta is list of length m with n x n matrices 

# b is list of length m with p x p matrices 

# x is n x p 

smacofIDIOSCAL <- 
function (w, 

delta, 

b, 

x, 

itmax = 100, 
eps = le-6, 
verbose = TRUE) { 
n <- nrow (x) 
p <- ncol (x) 
m <- length (w) 
itel <- 1 
sold <- Inf 
repeat { 
sa <- 0.0 
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v <- matrix (0, n * p, n * p) 
u <- matrix (0, n * p, n * p) 
for (k in l:m) { 

cmat <- tcrossprod (b [ [k]]) 

xmat <- x b [ [k] ] 

dmat <- as.matrix (dist (xmat)) 

vmat <- -w [[k]] 

diag(vmat) <- -rowSums(vmat) 

bmat <- -delta[[k]] * ifelse (dmat == 0, 0, 1 / dmat) 
diag(bmat) <- -rowSums(bmat) 
sa <- 

sa + sum (w[[k]] * (delta[[k]] - dmat) ~ 2) / 2.0 
v <- v + kronecker (cmat, vmat) 
u <- u + kronecker (cmat, bmat) 

} 

x <- matrix (ginv (v) 1*1 u 1*1 as.vector (x) , n, p) 

sb <- 0.0 

for (k in l:m) { 

xmat <- x 1*1 b [ [k]] 

dmat <- as.matrix (dist (xmat)) 

vmat <- -w [[k]] 

diag(vmat) <- -rowSums(vmat) 

hmat <- crossprod (x, vmat 1*1 x) 

bmat <- -delta[[k]] * ifelse (dmat == 0, 0, 1 / dmat) 
diag(bmat) <- -rowSums(bmat) 
gmat <- crossprod (x, bmat 1*1 x) 
sb <- 

sb + sum (w[[k]] * (delta[[k]] - dmat) ~ 2) / 2.0 
kmat <- ginv (hmat) 1*1 gmat 
for (j in l:p) { 

b[[k]][, j] <- kmat 1*1 b[[k]][, j] 

} 

} 

if (verbose) { 
cat ( 

"itel ", 

formate (itel, digits = 3, format = "d") , 

"sold ", 
formate ( 
sold, 

digits = 6, 
width = 10, 
format = "f" 

), 
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sa 


formate ( 

sa, 

digits = 6, 
width = 10, 
format = "f" 

), 

"sb ", 
formate ( 

sb, 

digits = 6, 
width = 10, 
format = "f" 

), 

"\n" 

) 

} 

if ((itel == itmax) |I (sold - sb < eps)) 
break 
sold <- sb 
itel <- itel + 1 

> 

return (list (x = x, b = b, s = sb)) 

} 

10.2 smacoflNDSCAL.R 

# w is list of length m with n x n matrices 

# delta is list of length m with n x n matrices 

# b is list of length m with p x p matrices 

# x is n x p 

smacofINDSCAL <- function (w, delta, b, x, itmax = 100, eps = le-6, verbose = TRUE) { 
n <- nrow (x) 
p <- ncol (x) 
m <- length (w) 
itel <- 1 
sold <- Inf 
repeat { 
sa <- 0.0 

v <- matrix (0, n * p, n * p) 
u <- matrix (0, n * p, n * p) 
for (k in l:m) { 

cmat <- tcrossprod (b [ [k]]) 
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xmat <- x 1*1 b [ [k] ] 

dmat <- as.matrix (dist (xmat)) 

vmat <- -w[[k]] 

diag(vmat) <- -rowSums(vmat) 

bmat <- -delta[[k]] * ifelse (dmat == 0, 0, 1 / dmat) 
diag(bmat) <- -rowSums(bmat) 
sa <- 

sa + sum (w[[k]j * (delta[[k]] - dmat) ~ 2) / 2.0 
v <- v + kronecker (cmat, vmat) 
u <- u + kronecker (cmat, bmat) 

> 

x <- matrix (ginv (v) 1*1 u 1*1 as.vector (x), n, p) 

sb <- 0.0 

for (k in l:m) { 

xmat <- x 1*1 b[[k]] 

dmat <- as.matrix (dist (xmat)) 

vmat <- -w[[k]] 

diag(vmat) <- -rowSums(vmat) 

hmat <- crossprod (x, vmat 1*1 x) 

bmat <- -delta[[k]] * ifelse (dmat == 0, 0, 1 / dmat) 
diag(bmat) <- -rowSums(bmat) 
gmat <- crossprod (x, bmat 1*1 x) 
sb <- 

sb + sum (w[[k]j * (delta[[k]] - dmat) ~ 2) / 2.0 
kmat <- ginv (hmat) °/ 0 *°/ 0 gmat 

b[[k]] <- diag (drop (kmat 1*1 diag (b[[k]]))) 

> 

if (verbose) { 
cat ( 

"itel ", 

formate (itel, digits = 3, format = "d"), 

"sold ", 
formate ( 
sold, 

digits = 6, 
width = 10, 
format = "f" 

), 

"sa ", 
formate ( 
sa, 

digits = 6, 
width = 10, 
format = "f" 
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), 

"sb ", 
formate ( 
sb, 

digits = 6, 
width = 10, 
format = "f" 

), 

"\n" 

) 

> 

if ((itel == itmax) I I (sold - sb < eps)) 
break 
sold <- sb 
itel <- itel + 1 

} 

return (list (x = x, b = b, s = sb)) 


10.3 smacofDerivatives.R 

smacofDerivativesX <- function (w, delta, b, x) 
n <- nrow (x) 
p <- ncol (x) 
m <- length (w) 
s <- 0.0 

v <- matrix (0, n * p, n * p) 

u <- matrix (0, n * p, n * p) 

r <- matrix (0, n * p, n * p) 

for (k in l:m) { 

cmat <- tcrossprod (b [ [k]]) 
xmat <- x y.*y, b [ [k] ] 
dmat <- as.matrix (dist (xmat)) 
for (i in 1: (n - 1)) { 
for (j in (i + 1):n) { 

ei <- ifelse (i == l:n, 1, 0) 

ej <- ifelse (j == l:n, 1, 0) 

amat <- outer (ei - ej, ei - ej) 
kmat <- kronecker (cmat, amat) 
kx <- drop (kmat °/ 0 *°/ 0 as .vector (x)) 
kxkx <- outer (kx, kx) 
s <- s + w[[k]][i, j] * (delta[[k]] [i, 
v <- v + w[ [k] ] [i, j] * kmat 


- dmat [i, j]) ~ 2 
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u <- 

u + w[[k]][i, j] * (delta[[k]] [i, j] / dmat[i, j]) * kmat 
r <- 

r + w[[k]][i, j] * (delta[[k]] [i, j] / (dmat [i, j] 3)) * kxkx 

> 

> 

> 

g <- drop ((v - u) 1*1 as .vector (x)) 
h <- v - u + r 

return (list (s = s, g = 0, h = h)) 

> 

smacofDerivativesTKDiag <- function (w, delta, b, x) { 
n <- nrow (x) 
p <- ncol (x) 
s <- 0.0 

v <- matrix (0, p, p) 

u <- matrix (0, p, p) 

r <- matrix (0, p, p) 

cmat <- tcrossprod (b) 

xmat <- x 1*1 b 

dmat <- as.matrix (dist (xmat)) 
for (i in 1: (n - 1) ) { 
for (j in (i + 1) :n) { 

ei <- ifelse (i == l:n, 1, 0) 
ej <- ifelse (j == l:n, 1, 0) 
amat <- outer (ei - ej, ei - ej) 
kmat <- crossprod(x, amat 1*1 x) 
kx <- drop (kmat 1*1 diag(b)) 
kxkx <- outer (kx, kx) 

s <- s + w[i, j] * (delta[i, j] - dmat [i, j]) ~ 2 
v<-v+w[i, j] * kmat 
u <- 

u + w[i, j] * (delta[i, j] / dmat[i, j]) * kmat 
r <- 

r + w[i, j] * (delta[i, j] / (dmat[i, j] “3)) * kxkx 

> 

} 

g <- drop ((v - u) 1*1 diag(b) ) 
h <- v - u + r 

return (list (s = s, g = g, h = h)) 

} 

smacofDerivativesTKFull <- function (w, delta, b, x) { 
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n <- nrow (x) 
p <- ncol (x) 
m <- length (w) 
s <- 0.0 

v <- matrix (0, n * p, n * p) 
u <- matrix (0, n * p, n * p) 
r <- matrix (0, n * p, n * p) 
cmat <- tcrossprod (b[[k]]) 
xmat <- x •/,*•/, b [ [k] ] 
dmat <- as.matrix (dist (xmat)) 
for (i in 1: (n - 1) ) { 
for (j in (i + 1) :n) { 

ei <- ifelse (i == l:n, 1, 0) 
ej <- ifelse (j == l:n, 1, 0) 
amat <- outer (ei - e j , ei - ej) 
kmat <- kronecker (cmat, amat) 
kx <- drop (kmat 1*1 diag(x)) 

s <- s + w[[k]][i, j] * (delta [ [k] ] [i, j] - dmat [i, j]) ~ 2 
v <- v + w[[k]] [i, j] * kmat 
u <- 

u + w[[k]][i, j] * (delta [ [k] ] [i, j] / dmat[i, j]) * kmat 
r <- 

r + w[[k]][i, j] * (delta [ [k] ] [i, j] / dmat[i, j]) * (outer (kx, kx) / (dmat [i, 

> 

} 

g <" drop ((v - u) 1*1 as.vector (x)) 
h <- v - u + r 

return (list (s = s, g = g, h = h)) 


10.4 smacofEllipses.R 

smacofEllipsesX <- function (x, h, f, eps, s = 1, t = 2) { 
n <- nrow (x) 
plot ( 
x, 

col = "RED", 
pch = 21, 
bg = "RED", 

xlab = paste ("dim ", formatC( 
s, digits = 1, format = "d" 

)), 

ylab = paste("dim ", formatC( 
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t, digits = 1, format = "d" 

)) 

) 

z <- 

cbind (sin (seq(-pi, pi, length = 100)), 

cos (seq(-pi, pi, length = 100))) * sqrt(2 * eps * f) 
for (i in l:n) { 

ii <- c ((s - 1) * n + i, (t - 1) * n + i) 
y <- x[i, c (s, t)] 
amat <- h[ii, ii] 
heig <- eigen (amat) 

zi <- z y o *°/ t diag (1 / sqrt (heig$values)) 
zi <- tcrossprod(zi , heig$vectors) 
zi <- zi + matrix (y, 100, 2, byrow = TRUE) 
lines (list(x = zi[, 1] , y = zi[, 2])) 

} 

} 

smacofEllipsesTDiag <- function (b, h, f, eps, s = 1, t = 2) { 
par (pty= "s") 
m <- nrow (b) 
ii <- c(s, t) 
plot ( 
b[, ii] , 
col = "RED", 
pch = 21, 
bg = "RED", 

xlab = paste ("dim ", formatC( 

s, digits = 1, format = "d" 

)), 

ylab = paste ("dim ", formatC( 

t, digits = 1, format = "d" 

)) 

) 

z <- 

cbind (sin (seq(-pi, pi, length = 100)), 

cos (seq(-pi, pi, length = 100))) * sqrt (2 * eps * f) 
for (k in l:m) { 
y <- b [k, ii] 
amat <- h[ii, ii, k] 
heig <- eigen (amat) 

zi <- z •/,*•/, diag (1 / sqrt (heig$values)) 

zi <- tcrossprod(zi , heig$vectors) 

zi <- zi + matrix (y, 100, 2, byrow = TRUE) 
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lines (list(x = zi[, 1] , y = zi[, 2])) 

> 

} 
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